##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--muti_cancer"), type = "character") ,
    make_option(c("--muti_pre"), type = "character") ,
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--type"), type = "character") ,
    make_option(c("--sample_class_file"), type = "character") ,
    make_option(c("--gene_list"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    muti_cancer <- "~/20220915_gastric_multiple/dna_combinePublic/maf/All_GGA.cancer.maf"
    muti_pre <- "~/20220915_gastric_multiple/dna_combinePublic/maf/All_GGA.precancer.maf"
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.tsv"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    gene_list <- "~/20220915_gastric_multiple/dna_combinePublic/mutsig_check/smg.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/GeneShare"
    type <- "IGC"
    sample_class_file <- "~/20220915_gastric_multiple/dna_combinePublic/images/GeneShare/Driver_Trunk.evolutionChoose.IGC.tsv"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene_list <- opt$gene_list
muti_cancer <- opt$muti_cancer
muti_pre <- opt$muti_pre
sample_info <- opt$sample_info
out_path <- opt$out_path
ccf_file <- opt$ccf_file
type <- opt$type
sample_class_file <- opt$sample_class_file

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_mutiCancer <- fread(muti_cancer  ,sep = "\t",  quote="" ,header = T)
dat_mutiPre <- fread(muti_pre  ,sep = "\t",  quote="" ,header = T)
dat_mutiNodule <- rbind(dat_mutiCancer , dat_mutiPre) 

dat_info <- data.frame(fread(sample_info))
dat_gene <- data.frame(fread(gene_list))
dat_ccf <- data.frame(fread(ccf_file))
dat_sampleInfo <- data.frame(fread(sample_class_file))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################
dat_mutiNodule <- subset(dat_mutiNodule , t_alt_count > 0)

dat_mutiNodule <- data.frame(Hugo_Symbol = dat_mutiNodule$Hugo_Symbol,
    Chromosome = dat_mutiNodule$Chromosome , Start_Position =  dat_mutiNodule$Start_Position , End_Position = dat_mutiNodule$End_Position ,
    Reference_Allele = dat_mutiNodule$Reference_Allele , Tumor_Seq_Allele2 = dat_mutiNodule$Tumor_Seq_Allele2 , 
    Variant_Classification = dat_mutiNodule$Variant_Classification,
    t_ref_count = dat_mutiNodule$t_ref_count , t_alt_count = dat_mutiNodule$t_alt_count ,
    Tumor_Sample_Barcode = dat_mutiNodule$Tumor_Sample_Barcode)

dat_mutiNodule$Location <- paste( dat_mutiNodule$Hugo_Symbol , dat_mutiNodule$Chromosome , dat_mutiNodule$Start_Position , 
    dat_mutiNodule$Reference_Allele , dat_mutiNodule$Tumor_Seq_Allele2 , sep=":"  )

dat_mutiNodule <- subset( dat_mutiNodule , Variant_Classification %in% Variant_Type)

###########################################################################################
## 注释突变的CCF
dat_mutiNodule$mergeCCF <- paste( dat_mutiNodule$Tumor_Sample_Barcode , dat_mutiNodule$Location , sep = ":" )

dat_ccf$mergeCCF <- paste( dat_ccf$Sample , dat_ccf$Hugo_Symbol , dat_ccf$Chr , dat_ccf$Start_Position , 
    dat_ccf$REF , dat_ccf$ALT , sep=":"  )

dat_mutiNodule <- merge( dat_mutiNodule , dat_ccf[,c("mergeCCF" , "minor_cn" , "total_cn" , "CCF_adj")] , by = "mergeCCF" )
dat_mutiNodule <- merge( dat_mutiNodule , dat_info[,c("Tumor" , "Class")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )

###########################################################################################
## 关注的克隆簇
if(type == "IGC"){
    dat_info <- subset( dat_info , Type == "IM + IGC" & Type != "IM + IGC + DGC" )
}else if(type == "DGC"){
    dat_info <- subset( dat_info , Type == "IM + DGC" & Type != "IM + IGC + DGC")
}else if(type == "All"){
    #dat_info <- subset( dat_info , Type != "IM + IGC + DGC" )
    dat_info <- dat_info
}


###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

## 判断基因多少比例出现在Trunk、Pre_Private、Inv_Private
result <- c()
result_plot <- c()
gene_list <- dat_gene$Gene_Symbol

for( geneN in gene_list ){

    tmp_info <- subset( dat_sampleInfo , gene == geneN )

    shareSample <- unlist(strsplit(c(tmp_info$trunk_driver_sample , tmp_info$share_driver_sample) , ","))
    invSample <- unlist(strsplit(c(tmp_info$private_driver_inv_sample) , ","))

    shareLOH_num <- 0
    for( Sample in shareSample ){
        tumors <- subset( dat_info , ID == Sample )$Tumor
        tmp <- subset( dat_mutiNodule , Tumor_Sample_Barcode %in% tumors & Variant_Classification %in% Variant_Type & Hugo_Symbol == geneN )
        if( nrow(subset( tmp , minor_cn == 0 )) > 0 ){
            shareLOH_num <- shareLOH_num + 1 
        }
    }
    shareNonLOH_num <- length(shareSample) - shareLOH_num

    privateLOH_num <- 0
    for( Sample in invSample ){
        tumors <- subset( dat_info , ID == Sample )$Tumor
        tmp <- subset( dat_mutiNodule , Tumor_Sample_Barcode %in% tumors & Variant_Classification %in% Variant_Type & Hugo_Symbol == geneN )
        if( nrow(subset( tmp , minor_cn == 0 )) > 0 ){
            privateLOH_num <- privateLOH_num + 1 
        }
    }
    privateNonLOH_num <- length(invSample) - privateLOH_num

    dat_matrix <- matrix(c( shareLOH_num , shareNonLOH_num , privateLOH_num , privateNonLOH_num ), ncol = 2)
    p <- fisher.test(dat_matrix)$p.value
    or <- fisher.test(dat_matrix)$estimate

    if( p < 0.001 ){
        p_text <- trans(p)
    }else{
        p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
    }
    
    tmp_use <- data.frame( gene = geneN , shareLOH_num = shareLOH_num , shareNonLOH_num = shareNonLOH_num , 
        privateLOH_num = privateLOH_num , privateNonLOH_num = privateNonLOH_num , p = p , OR = or ) 

    result <- rbind( result , tmp_use )

    tmp_use <- data.frame( gene = geneN , 
        type = c("Share" , "Share" , "Private" , "Private") ,
        Num = c(shareLOH_num , shareNonLOH_num , privateLOH_num , privateNonLOH_num) ,
        Ratio = c( 
            shareLOH_num/(shareLOH_num+shareNonLOH_num) , shareNonLOH_num/(shareLOH_num+shareNonLOH_num) ,
            privateLOH_num/(privateLOH_num+privateNonLOH_num) , privateNonLOH_num/(privateLOH_num+privateNonLOH_num) 
            ) ,
        loh = c("LOH" , "Non-LOH" , "LOH" , "Non-LOH")
     , p = p , OR = or ) 
    tmp_use$p_text <- ""
    tmp_use$p_text[1] <- p_text
    tmp_use$OR <- ""
    tmp_use$OR[1] <- or

    result_plot <- rbind(result_plot , tmp_use)
}


out_name <- paste0(out_path , "/Driver_Trunk.",type,".LOH.tsv")
write.table( result , out_name , row.names = F , sep = "\t" , quote = F )

###########################################################################################
## 基因分布堆叠图
dat <- result_plot
dat$Ratio[is.na(dat$Ratio)] <- 0
dat$value_text <- paste0( round(dat$Ratio , 2) * 100 , "%")

plot <- ggplot( data = dat , aes( x = type , y = Ratio , fill = loh ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion")+
    facet_grid(.~gene) +
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='bottom',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold',angle = 45, vjust = 1, hjust=1) ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/Driver_Trunk.",type,".LOH.pdf")
ggsave( out_name , plot , width = 15 , height = 5 )

